An efficient, multiple range random walk algorithm to calculate the density of states 
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We present a new Monte Carlo algorithm that produces results of high accuracy with reduced 
simulational effort. Independent random walks are performed (concurrently or serially) in different, 
restricted ranges of energy , and the resultant density of states is modified continuously to produce 
locally flat histograms. This method permits us to directly access the free energy and entropy, is 
independent of temperature, and is efficient for the study of both 1st order and 2nd order phase 
transitions. It should also be useful for the study of complex systems with a rough energy landscape. 
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Computer simulation has become an essential tool in 
condensed matter physics JlJ , particularly for the study of 
phase transitions and critical phenomena. The workhorse 
for the past half-century has been the Metropolis impor- 
tance sampling algorithm, but more recently new, effi- 
cient algorithms have begun to play a role in allowing 
simulation to achieve the resolution which is needed to 
accurately locate and characterize phase transitions. For 
example, cluster flip algorithms, beginning with the sem- 
inal work of Swendsen and Wang [Q, have been used 
to reduce critical slowing down near 2nd order transi- 
tions. Similarly, the multicanonical ensemble method ||] 
was introduced to overcome the tunneling barrier be- 
tween coexisting phases at 1st order transitions, and this 
approach also has utility for systems with a rough en- 
ergy landscape ||-^|- In both situations, histogram re- 
weighting techniques R| can be applied in the analysis to 
increase the amount of information that can be gleaned 
from simulational data, but the applicability of reweight- 
ing is severely limited in large systems by the statistical 
quality of the "wings" of the histogram. This latter effect 
is quite important in systems with competing interactions 
for which short range order effects might occur over very 
broad temperature ranges or even give rise to frustration 
that produces a very complicated energy landscape and 
limit the efficiency of other methods. 

In this paper, we introduce a new, general, efficient 
Monte Carlo algorithm that offers substantial advantages 
over existing approaches. Unlike conventional Monte 
Carlo methods that directly generate a canonical dis- 
tribution at a given temperature g(E)e~ E ^ KliT , our ap- 
proach is to estimate the density of states g{E) accurately 
via a random walk which produces a flat histogram in en- 
ergy space. The method can be further enhanced by per- 
forming multiple random walks, each for a different range 
of energy, cither serially or in parallel fashion. The resul- 
tant pieces of the density of states can be joined together 
and used to produce canonical averages for the calcula- 
tion of thermodynamic quantities at essentially any tem- 
perature. We will apply our algorithm to the 2-dim ten 
state Potts model and Ising model which have 1st- and 
2nd-order phase transitions, respectively, to demonstrate 



the efficiency and accuracy of the method. 

Our algorithm is based on the observation that if we 
perform a random walk in energy space with a probabil- 
ity proportional to the reciprocal of the density of states 
g ^~ E j , then a flat histogram is generated for the energy 
distribution. This is accomplished by modifying the es- 
timated density of states in a systematic way to produce 
a "flat" histogram over the allowed range of energy and 
simultaneously making the density of states converge to 
the true value. At the very beginning of the random 
walk, the density of states is a priori unknown, so we 
simply set all densities of states g(E) for all energies E 
to g(E) = 1. Then we begin our random walk in en- 
ergy space by flipping spins randomly. In general, if E\ 
and Ei are energies before and after a spin is flipped, 
the transition probability from energy level E± to E% is 
simply: 



p(Ei — > E2) = min( 



g(Ei) 
g{E 2 y 
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This is also the probability to flip the spin. Each time 
an energy level E is visited, we update the corresponding 
density of states by multiplying the existing value by a 
modification factor / > 1, i.e. g(E) — > g(E) * f. The ini- 
tial modification factor can be as large as / = fo = e 1 ~ 
2.71828... which allows us to reach all possible energy lev- 
els very quickly, even for large systems. We keep walking 
randomly in energy space and modifying the density of 
states until the accumulated histogram H(E) is "flat" . 
At this point, the density of states converges to the true 
value with an accuracy proportion to ln(/) . We then 
reduce the modification factor to a finer one according 
to some recipe like f± — vJo (any function that mono- 
tonically decreases to 1 will do) and reset the histogram 
H(E) = 0. Then we begin the next level random walk 
with a finer modification factor / = /1 , continuing until 
the histogram is again "flat" after which we stop and re- 
duce the modification factor as before, i.e. = \ffi- 
We stop the simulation process when the modification 
factor is smaller than some predefined final value (such 
as / fina i = exp(10~ 8 ) ~ 1.00000001). It is very clear 
that the modification factor / in our random walk acts 



1 



as a control parameter for the accuracy of the density 
of states during the simulation and also determines how 
many MC sweeps are necessary for the whole simulation. 
It is impossible to obtain a perfectly flat histogram and 
the phrase "flat histogram" in this paper means that his- 
togram H(E) for all possible E is not less than 80% of the 
average histogram (H{E)). Since the density of states is 
modified every time the state is visited, we only obtain a 
relative density of states at the end of the simulation. To 
calculate the absolute values, we use the condition that 
the number of ground states for the Ising model is 2 (all 
spins are up or down) to re-scale the density of states; 
and if multiple walks are performed within different en- 
ergy ranges, they must be matched up at the boundaries 
in energy. 

Because of the exponential growth of the density of 
states in energy space, it is not efficient to simply up- 
date the density of states until enough histogram en- 
tries are accumulated. All methods based on the ac- 
cumulation of entries, such as the histogram method 

0, Lee's version of the multicanonical method (entropic 
sampling) ||, the broad histogram method Q and the 
flat histogram method P]lO|] have the problem of scala- 
bility for large systems. These methods suffer from sys- 
tematic errors and substantial deviations which increase 
rapidly for large system size. The algorithm proposed in 
this paper is of both high efficiency and accuracy over 
wide ranges of temperature for sizes that are beyond 
those that are tractable by other approaches. 

We should point out here that during the random walk 
(especially for the early stage of iteration) , the algorithm 
does not exactly satisfy the detailed balance condition, 
since the density of states is modified constantly during 
the random walk in energy space; however, after many 
iterations, the density of states converges to the true 
value very quickly as the modification factor approaches 

1. From eq. (1), we have: 
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g(Ei) 



p{E x -» E 2 ) 



1 



g(E 2 ) 



P (E 2 -» Ex) 



(2) 



where 



i s t ne probability at the energy level E\ and 
•+ E 2 ) is the transition probability from E\ to E 2 
We can thus conclude that the 



p(Ei 

for the random walk 
detailed balance condition is satisfied to within the accu- 
racy proportion to ln(/). 

The convergence and accuracy of our algorithm may 
be tested for a system with a 2nd order transition, the 
L x L Ising square lattice with nearest neighbor coupling 
which is generally perceived as an ideal benchmark for 
new theories and simulation algorithms We 
simulated both small lattices for which exact results are 
available as well as L = 256 for which exact enumeration 
is impossible. In Fig. 1, the densities of states estimated 
by our algorithm are shown along with the exact results 
obtained by the method proposed by Beale Q . We only 



show the density for systems up to L — 50 which is the 
maximum size we can calculate with the Mathematica 
program used in the reference |l3| . Since no difference is 
visible, we show the relative error e(\og(g(E))), which is 
defined as e(X) = \(X sim - X cxact ) / X cxact \ for a general 
quantity X in this paper. With our algorithm we ob- 
tain an average error as small as 0.035 % on the 32 x 32 
lattice with 7 x 10 5 sweeps. It is possible to estimate 
the density of states for small systems with the broad 
histogram method S. Recent broad histogram simula- 
tional data H] for the 2D Ising model on a 32 x 32 lattice 
with 10 6 MC sweeps yielded an average deviation of the 
microcanonical entropy from about 0.08 % from the exact 
solution (l3|. 
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FIG. 1. Comparison of the density of states obtained 
by our algorithm for 2D Ising model and the exact results 
calculated by the method in reference [13]. Relative errors 
(e(\og(g(E))) are shown in the inset. 

With the Monte Carlo algorithm proposed in this pa- 
per, we can estimate the density of states efficiently even 
for large systems. Because of the symmetry of the density 
of states for Ising model g{E) = g(—E), we only need to 
estimate the density of states in the region E/N £ [—2, 0] , 
where N is total lattice sites. To speed up the simulation 
for L = 256, we perform 15 independent random walks, 
each for a different region of energy from E/N = —2 
to E/N = 0.2 using f &na[ = exp(10~ 8 ). To reduce the 
"boundary effect", random walks over adjacent energy 
regions overlap by AE/N — 0.06. The density of states 
for E/N 6 [—2,0.2] is obtained by joining 15 densities 
of states from random walks on different energy regions 
using a total simulational effort of only 6.1 x 10 6 MC 
sweeps. 

One advantage of our algorithm is that we can readily 
calculate the Gibbs free energy and the entropy, quan- 
tities which are not directly available in conventional 
Monte Carlo simulations. With the density of states, 
the Gibbs free energy can be calculated by 
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F{T) = -k B Tln(Z) 



-k B Tln(Y,9(E)e- 

E 



■PE 



)• (3) 



Although it is impossible to calculate the exact density of 
states of Ising model on a lattice as large as L — 256 with 
the method proposed by Beale p3| , the free energy and 
specific heat were calculated exactly by Ferdinand and 
Fisher on finite-size lattices. In Fig. 2, we compare 
simulational data and exact solutions for the Gibbs free 
energy as a function of temperature. The agreement is 
excellent and a more stringent test of the accuracy shows 
that the relative error e(F) is smaller than 0.0008% for 
temperature region T e [0, 8]. 
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Our simulational data on the finite-size lattice are com- 
pared with the exact solution obtained by Ferdinand 
and Fisher jH in Fig. 3. A str ingent test of the ac- 
curacy is provided by the inset which shows the rela- 
tive error e(C). The average error over the entire range 
T € [0.4,8] only used a total of 6.1 x 10 6 MC sweeps 
is 0.39%. The relative errors are not bigger than 4.5% 
even with fine scale near T c . Recently, Wang, Tay and 
Swendsen estimated the specific heat of the same 
model on a 64 x 64 lattice by the transition matrix Monte 
Carlo re- weighting method JItJ , and for a simulation with 
2.5 x 10 7 MC sweeps, the maximum error in temperature 
region T S [0, 8] was about 1%. When we apply our al- 
gorithm to the same model on the 64 x 64 lattice, with 
a final modification factor of 1.000000001 and a total of 
2 x 10 7 MC sweeps on single processor, the errors of the 
specific heat are reduced below 0.7% for all temperature 
p6[ . The relatively large errors at low temperature reflect 
the small values for the specific heat at low temperature. 
The errors in specific heat estimated from the density of 
states with broad histogram method are obviously visible 
even for systems as small as 32 x 32 || whereas with our 
method, such differences are invisible even for a system 
as large as 256 x 256. 
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FIG. 2. Comparison of the Gibbs free energy per lattice 
site calculated directly from the density of states from our 
simulation for the L = 256 Ising model and the exact solutions 
from reference [15]. The relative errors s(F) are shown in the 
inset. The density of states was obtained by random walks 
with only 6.1 x 10 6 MC sweeps totally. 

The entropy is a very important thermodynamic quan- 
tity that cannot be calculated directly in conventional 
Monte Carlo simulations. It can be estimated by in- 
tegrating over other thermodynamic quantities, such 
as specific heat, but such calculations are not so reli- 
able since the specific heat itself is not so easy to es- 
timate accurately. With an accurate density of states 
estimated by our method, the entropy can be calcu- 
lated easily by S{T) = u ( T )- p( ~ T ) wne re U(T) = 
(E) T = Eg(E)e~ fJE /J2 g(E) e - f}E is the internal en- 

E E 

ergy. According to our calculation, the errors for L = 
256 are smaller than 1.2% in all temperature region 
T 6 [0,8]. |l6| Very recently, with the flat histogram 
method || and the broad histogram method |^| , the en- 
tropy was estimated with 10 7 MC sweeps for the same 
model on 32 x 32 lattice; however, the errors in reference 
JTo| are even much bigger than our errors for 256 x 256! 

A more stringent test of the accuracy of the density 
of states is calculation of the specific heat defined by the 
fluctuation expression: 




FIG. 3. Specific heat for the 2-dim Ising model on a 
256 x 256 lattice in a wide temperature region. The relative 
error e(C) are shown in the inset in the figure. 

With our algorithm, we not only dramatically reduce 
the computational effort by avoiding multiple simulations 
for different temperatures close to the transition, but also 
overcome the slow kinetics at low temperature or near T c 
for both first-order and second-order phase transitions 
since the random walk does not depend on the tempera- 
ture. To show how our simulation method overcomes the 
tunneling barrier between order and disorder phases at 
a lst-order phase transition, we perform random walks 



3 



to calculate the density of states for the 2D ten state 
Potts model |Q with nearest neighbor interactions on 
square lattices of size 60 < L < 200. In Fig. 4, we show 
the canonical distributions at the temperatures at which 
the peaks are of equal height. Because of the double peak 
structure of strongly lst-order phase transitions |Q, con- 
ventional Monte Carlo simulations are not efficient since 
it takes an extremely long time to tunnel from one peak 
to the other. Considering the valley which we find for 
L = 200 is as deep as 9 x 10~ 10 , it is impossible for 
conventional Monte Carlo algorithms to overcome such a 
tunneling barrier with available computational resources. 



should be extremely useful for the study of the complex 
systems such as spin glass models jSJ and protein folding 
problems |6| where the energy landscape is very rough 
and where it is already known that there are problems 
with other optimization algorithms. 

We would like to thank C. K. Hu, N. Hatano, P. D. 
Beale, S. P. Lewis and H-B Schuttler for helpful com- 
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ported by the National Science Foundation under Grant 
No. DMR-9727714. 
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FIG. 4. The canonical distributions for the 2-dim ten state 
Potts model on L x L lattice at the transition temperature 
P{E,T C ) = g(E)e~ E/KBT °. For L = 150 and 200, multiple 
random walks were performed in different energy regions with 
locally flat histograms. The distributions at peaks are nor- 
malized to 1. The transition temperature T C (L) is 0.701243 
for L = 200; T c (oo) = 0.701232. ...(exact solution). In the 
inset, we show the overall histograms of the random walks 
for L = 100 and 200. 3.1 x 10 7 visits per energy level were 
used for L = 100 with a single random walk. With multiple 
random walks, the density of states for L — 200 was obtained 
with only 9.8 x 10 6 visits per energy level. 

All thermodynamic quantities we discussed so far are 
directly related to energy. It is also possible to cal- 
culate any quantities which may not directly relate to 
energy. |l(| As an example, the order-parameter for 
the 2D Ising model can be calculated by |Af(T)| = 
J2 \M{E)\g(E)e-P E /J2 g{E)e~ t}E where M(E) is the av- 

E E 

erage value of the order-parameter at energy level E dur- 
ing the random walk. The random walk is not restricted 
to energy space, and our algorithm can be applied to any 
other parameter space. To apply our algorithm to a new 
system, the only thing we need to know is the Hamilto- 
nian. The algorithm can be optimized to estimate the 
relevant density of states to the property and tempera- 
ture range of interest. This new Monte Carlo algorithm 



[6 
[7; 
[8 

[9 
[10 

[11 
[12 

[13 
[14 

[IS 

[16 
[17 

[18 
[19 



D. P. Landau and K. Binder, A Guide to Monte Carlo 
Methods in Statistical Physics, (Cambridge U. Press, 
Cambridge, 2000). 

R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 

(1987), U. Wolff, Phys. Rev. Lett. 62, 361 (1989). 

B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 

(1992), Phys. Lett. B 267, 249 (1991), J. Lee, Phys. Rev. 

Lett. 71 , 211 (1993). B. A. Berg J. Stat. Phys. 82, 323 

(1996), B. A. Berg, Nucl. Phys. B, 63, 982 (1998). 

W. Janke and S. Kappler, Phys. Rev. Lett. 74, 212 

(1995). 

B. A. Berg and T. Celik, Phys. Rev. Lett. 69, 2292 
(1992), B. A. Berg and W. Janke, Phys. Rev. Lett. 80, 
4771 (1998), N. Hatano and J. E. Gubernatis in Com- 
puter Simulation Studies in Condensed Matter Physics 
XII, D. P. Landau, S. P. Lewis and H.-B. Schuttler (eds) 
(Springer, Berlin, Heidelberg, 2000). 
N. S. Alves and U. Hansmann, Phys. Rev. Lett. 84, 1836 
(2000). 

A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 
61, 2635 (1988), 63, 1195 (1989). 

P. M. C. de Oliveira, T. J. P. Penna and H. J. Herrmann, 
Eur. Phys. J. B. 1, 205 (1998), P. M. C. de Oliveira, Eur. 
Phys. J. B. 6, 111 (1998). 

Jian-Sheng Wang, Eur. Phys. J. B. 8, 287 (1998), Cond- 
matt-9909177. 

Jian-Sheng Wang and Lik Wee Lee, Comput. Phys. Com- 
mun. 127, 131, (2000). 

D. P. Landau, Phys Rev. B 13, 2997 (1976). 

Jian-Sheng Wang, T. K. Tay and R. H. Swendsen, Phys. 

Rev. Lett. 82, 476 (1999). 

P. D. Beale, Phys. Rev. Lett. 76, 78 (1996). 

A. R. Lima, P. M. C. de Oliveira and T. J. P. Penna, J. 

Stat. Phys. 99, 691, (2000). 

A. E. Ferdinand and M. E. Fisher, Phys. Rev. 185, 832 
(1969). 

Fugao Wang and D. P. Landau, unpublished data. 

R. H. Swendsen, J. S. Wang, S. T Li, B. Diggs, C. Gen- 

ovese and J. B. Kadane, Cond-matt-9908461. 

F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982). 

M. S. S. Challa, D. P. Landau and K. Binder Phys. Rev. 

B 34, 1841, (1986). 



4 



